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Abstract. We analyze the Bethe ansatz equations describing the complete spectrum 
of the transition matrix of the partially asymmetric exclusion process on a finite lattice 
and with the most general open boundary conditions. We extend results obtained 
recently for totally asymmetric diffusion [J. de Gier and F.H.L. Essler, J. Stat. Mech. 
P12011 (2006)] to the case of partial asymmetry. We determine the finite-size scaling of 
the spectral gap, which characterizes the approach to stationarity at late times, in the 
low and high density regimes and on the coexistence line. We observe boundary induced 
crossovers and discuss possible interpretations of our results in terms of effective domain 
wall theories. 
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1. Introduction 

The partially asymmetric simple exclusion process (PASEP) [1,2] is one of the most 
thoroughly studied models of non-equilibrium statistical mechanics [3-6]. It is a 
microscopic model of a driven system [7] describing the asymmetric diffusion of hard- 
core particles along a one-dimensional chain with L sites. At late times the PASEP 
exhibits a relaxation towards a non-equilibrium stationary state. In the presence of 
two boundaries at which particles are injected and extracted with given rates, the 
bulk behaviour at stationarity is strongly dependent on the injection and extraction 
rates. The corresponding phase diagram as well as various physical quantities have 
been determined by exact methods [3,4,6,8-13]. 

Given the behaviour in the stationary state an obvious question is how the system 
relaxes to this state at late times. For the PASEP on a ring, where particle number is 
conserved, such results were obtained by means of Bethe's ansatz some time ago [14-16] . 
More recently there has been considerable progress in analyzing the dynamics in the limit 
of totally asymmetric exclusion and on an infinite lattice, see e.g. [17-22], where random 
matrix techniques can be used. 

For the finite system with open boundaries there have been several studies of 
dynamical properties by means of numerical, phenomenological and renormalization 
group methods [23-28]. For the case of symmetric diffusion a Bethe ansatz solution was 
constructed in [29]. Recently, we have applied Bethe's ansatz to the PASEP with open 
boundaries [30,31]. It is well-known that the PASEP can be mapped onto the spin-1/2 
anisotropic Heisenberg chain with general (integrable) open boundary conditions [10,11]. 
By building on recent progress in applying Bethe's ansatz to the latter problem [32-37] 
we determined the Bethe ansatz equations for the PASEP with the most general open 
boundary conditions. By analyzing these equations we derived the finite size scaling 
behaviour of the spectrum of low-lying excited states for the cases of symmetric and 
totally asymmetric diffusion. Upon varying the boundary rates, we observed crossovers 
in massive regions, with dynamic exponents z — 0, and between massive and scaling 
regions with diffusive [z = 2) and KPZ (z — 3/2) behaviour. 

In the present work we extend these results to the case of partially asymmetric 
diffusion, where the analysis of the spectrum is significantly more involved. 
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Figure 1. Transition rates for the partially asymmetric exclusion process. 



We now turn to a description of the dynamical rules defining the PASEP on a one 
dimensional lattice with L sites. At any given time t each site is either occupied by a 
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particle or empty. The system is then updated as follows. In the bulk of the system, 
i.e. sites % — 2, . . . , L — 1, a particle attempts to hop one site to the right with rate p 
and one site to the left with rate q. The hop is prohibited if the neighbouring site is 
occupied. On the first and last sites these rules are modified. If site % — 1 is empty, a 
particle may enter the system with rate a. If on the other hand site 1 is occupied by a 
particle, the latter will leave the system with rate 7. Similarly, at % — L particles are 
injected and extracted with rates 5 and j3 respectively. 

It is customary to associate a Boolean variable with every site, indicating whether 
a particle is present (7$ = 1) or not (7$ = 0) at site %. Let |0) and |1) denote the 
standard basis vectors in C 2 . A state of the system at time t is then characterized by 
the probability distribution 

\P(t)) = J2P(r\t)\r), (1.1) 

where 

L 

|T) = |n,...,r L ) = (g)|r i ). (1.2) 

i=i 

The time evolution of \P(t)) is governed by the aforementioned rules, which gives rise 
to the master equation 
d 

dt 1 

where the PASEP transition matrix M consists of two-body interactions only and is 
given by 

L-l 

M = I [k ~ l) ® M <g> J^-*" 1 ) + mi ® I^'^ + I {L - 1] ® m L . (1.4) 
fe=i 

Here J^-* is the identity matrix on the k-fold tensor product of C 2 and M : C 2 ® C 2 — > 
C 2 ® C 2 is given by 



\P(t)) = M\P(t)), (1.3) 
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The terms involving mi and m^ describe injection (extraction) of particles with rates a 
and 5 (7 and /3) at sites 1 and L respectively. Their explicit forms are 

.;)■ -=("/ -")■ (L6) 

The transition matrix M has a unique stationary state corresponding to the 
eigenvalue zero. For positive rates, all other eigenvalues of M have non-positive real 
parts. The late time behaviour of the PASEP is dominated by the eigenstates of M 
with the largest real parts of the corresponding eigenvalues, as follows from the following 
argument. The average of an observable X is given by (see e.g. Ref. [38]) 

(X)(t) = (0\Xe Mt \P ) (1.7) 
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Here Pq is an initial state and (0| is the left eigenstate of M with eigenvalue 0. This 
may be written in the spectral representation with respect to the eigenstates of M 

(X){t) = ^(0\X\n)e £ » t a n , (1.8) 

n 

where \Pq) = a n\ n ) and M\n) = S n \n). In the limit t — > oo only the stationary state 
survives and we have 

lim(X)(t) = (0|X|0)a . (1.9) 

t — >oo 

At very late times we have 

(X)(t) « (0|X|0)a + (0|X|l)a ie ^, (1.10) 

where E\ is the eigenvalue of M with the largest real part that contributes to the spectral 
decomposition of the initial state \Po)- Hence the relaxation rate —E\ determines the 
approach to the stationary state at asymptotically late times. In the next sections we 
determine the eigenvalue of M with the largest non-zero real part using Bethe's ansatz. 
The latter reduces the problem of determining the spectrum of M to solving a system 
of coupled polynomial equations of degree 3L — 1. Using these equations, the spectrum 
of M can be studied numerically for very large L, and, as we will show, analytic results 
can be obtained in the limit L — > oo. 



2. Bethe ansatz equations 

In [30,31] it was shown that the PASEP transition matrix M can be diagonalised using 
the Bethe Ansatz. In [31] the Bethe equations were analyzed in some detail for the case 
q = 7 = 5 = corresponding to totally asymmetric diffusion (TASEP). The resulting 
TASEP dynamical phase diagram displays interesting crossovers within the low and high 
density phases, at which the transition matrix eigenvalue corresponding to the slowest 
relaxation mode changes non-analytically. 

This work is concerned with the generalization of some of these results to the 
PASEP case with q, 7, 5 7^ 0. Before turning to the technical details of our analysis we 
present a summary of our main results. Throughout this work we set without loss of 
generality q < p = 1. For simplicity we only consider L even, as for odd L the details 
will be somewhat different. 

As was shown in [30,31], all eigenvalues £ of M can be expressed in terms of the 
roots Zj of a set of L — 1 non-linear algebraic equations as [|] 

g = -fi.-X L (q ~f Z \y (2-D 

where 

£ = a + f3 + 1 + 5. (2.2) 

I We rescale the roots zj in eqns (3.1) - (3.3) of [31] by Q = Jq. 
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The complex roots Zj satisfy the Bethe ansatz equations 



and 



qzj 



1 — Zi 



n 



qzj - Zj q^ZjZi 
z j ~ Q z i z j z i ~ 



Here K(z) = K(z,a,j)K(z, (3,8), where 

(z + K+ y )(z + K~ 7 ) 



(qK+„z + l){qK~ a z + 1) 



Q i7 



2a 

1 — q — a + 7. 



3 = 1. ..L 



In order to ease notations we will use the following abbreviations, 

a = K a,~f> b = K /3,<5> c = K a,-yi d = K /3,<5' 

The constant Sq is expressed in our new notations as 



£0 = (1 - q) 



1 — ac 



+ 



1 - bd 



l + a)(l + c) (l + 6)(l + d) 



(2.3) 

(2.4) 

(2.5) 
(2.6) 

(2.7) 
(2.8) 



3. Lowest excitation of the PASEP in the "forward-bias regime": summary 
of main results 



By analysing the set of equations (2.3) for large, finite L we have determined the 



eigenvalue of the transition matrix with the largest non-zero real part. From this "lowest 
excited state energy" we can infer properties of the relaxation towards the stationary 
state at asymptotically late times. In the present work we have restricted our analysis 
to the regime of small values of the parameters c and d. This corresponds loosely 
to a "forward-bias regime" in which particles diffuse predominantly from left to right 
and particle injection and extraction occurs mainly at sites 1 and L respectively. The 
restrictions on the permitted values of c and d are discussed in more detail in section [6j 



3.1. Stationary state phase diagram 

The phase diagram of the PASEP at stationarity was found by Sandow [10] and is 
depicted in Figure |2j We note that the phases depend only on the parameters a and b 



defined in (2.7) rather than p, q, a, f3, 7, 5 separately. 



3.2. Dynamical phase diagram 

The dynamical phase diagram for the PASEP resulting from our analysis in the regime 
q < p = 1 is shown in Figure [3j It exhibits the same subdivision into low and high 
density phases (a > 1 and b > 1), the coexistence line (a = b > 1) and the maximum 
current phase (a, b < 1) found from the analysis of the current in the stationary 



Slowest relaxation mode of the PASEP with open boundaries 6 

b 




MC 



low 

density 



Figure 2. Stationary state phase diagram of the PASEP. The high and low density 
phases are separated by the coexistence line (CL). The maximum current phase (MC) 
occurs at small values of the parameters a and b defined in (2.7 1. 



state [10]. However, the finite-size scaling of the lowest excited state energy of the 
transition matrix suggests the sub-division of both low and high-density phases into 
four regions respectively. These regions are characterized by different functional forms 
of the relaxation rates at asymptotically late times. We will describe the four regions 
in the low density phase (a > 1 and b < a). The results in the high density phase 
are obtained by exchanging a «-> b. We note that the results for the relaxation rates 
presented below are valid in the limit L — > oo at fixed q < 1. In particular the limit 
of symmetric exclusion q — 1 cannot be obtained by taking the limit q — > 1 in the 
expression presented below. 

(i) Region I: 

This region is defined by 

b > q 2 a for a > q~ 3/2 and b > a~ 1/3 for 1 < a < q~ 3/2 . (3.1) 
The eigenvalue of the lowest excitation is given by (q < 1 fixed) 



E x = - (1 - q) 



1 1 2z c 

+ 



1 + a 1 + 6 1 - z c 



+ (1 J t J 2 ,, +0(L-\ (3.2) 



i 2 (*„" 



where 



1 



ab 



(3.3) 



(ii) Coexistence Line: 

The coexistence line is defined by a = b > 1 and separates the low and high 



density phases. We find that the leading term in (3.2) vanishes and that the lowest 
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Si 



Q 



7T 



L 2 



(a 1 — a' 



+ 0{L 



-3n 



(3.4) 



The inverse proportionality of the eigenvalue (3.4) to the square of the system size 
implies a dynamic exponent z — 2, which in turn suggests that the relaxation at 
late times is governed by diffusive behaviour, 
(hi) Region II: 

This region is defined by 



b < cT 1/3 for 1 < a < q~ 3/2 



(3.5) 



The eigenvalue of the lowest excitation is now independent of b 
S x =-(\-q) 



1 2z c + 1 
+ 



+ 



1-q 



1 + a 
4tt 2 



+ 0{L~ 



where 



-a 



-1/3 



(3.6) 



(3.7) 



We note that the leading terms of (3.2) and (3.6) coincide along the boundary 
b = a -1 / 3 separating the two regimes, but the terms of order L~ 2 exhibit a 
discontinuity. This suggests a crossing of levels and an associated change in the 
detailed nature of the corresponding relaxational mode. 

(iv) Region III: 

This region is defined by 

q 1 / 2 <b< q 2 a for a > q~ 3 / 2 . (3.8) 

Up to terms of order 0(L~ 3 ), the eigenvalue of the lowest excitation in this region 
is given by 



Si 



q) 



1 



+ 



qa 



+ 



2£c-l 



1 + a 1 + b l + qa 1 



where now 



-q/b. 



(3.9) 



(3.10) 



We note that the leading terms of (3.2) and (3.9) coincide along the boundary 
b = q 2 a separating Regions I and III. However, throughout Region III there is 
no contribution of order 0(L~ 2 ) to the transition matrix eigenvalue of the lowest 
excited state. 

(v) Region IV: 

This final region is defined by 

a>q- 3/2 } b<q 1/2 . (3.11) 
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The eigenvalue of the lowest excitation in this region is given by 

1 qa 2z r 



Si 



+ 



1-q 



1 + a 

It 2 



where now 



-q 



L 2 (z^-z c ) 



1/2 



+ i 

1 + qa 1 



(3.12) 



(3.13) 



We note that the leading terms of (3.12) and (3.9) match along the boundary 
between regions IV and III. The same holds for the leading terms of (3.12) and 
(3.6) along the boundary between regions IV and II. On the other hand, there is a 
discontinuity in the 0(L~ 2 ) contributions in both cases. 




Figure 3. Dynamical phase diagram of the PASEP in the forward bias regime 
determined by the lowest excitation E\. The horizontal axes are the boundary 
parameters a and b (2.7) and the vertical axis is the lowest relaxation rate E\. The 
latter goes to zero for large systems on the coexistence line (CL) and in the maximum 
current phase (MC). The curves and lines correspond to various crossovers in the low 
and high density phase, across which E\ changes non-analytically. See the main text 
for a detailed explanation. 



3.3. Modified domain wall theory 

It was shown in [26, 27] that the diffusive relaxation towards the stationary state found 
on the coexistence line as well as in the low and high density phases can be understood 
in terms of an effective domain wall theory (DWT). In this approach the excited states 
driving the relaxational dynamics are modelled as domain walls between low and high 
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density regions. They carry out a random walk with right and left hopping rates given 
by 



;i-g) 



p 



(3.14) 



p^ - p 

Here, p~ = 1/(1 + a) and p + = b/(l + b) are the stationary bulk densities in the low and 
high-density phases respectively. The domain walls are assumed to be reflected from 
both boundaries. Interestingly, domain wall theory gives the exact stationary state along 
the curve ab = q^ 1 in parameter space [45]. It is furthermore possible to construct an 
entire family of exact domain wall solutions of the master equation [45] [§] The leading 
relaxation rate calculated from these exact domain wall solutions agrees with our results 



(3.2) and (3.4) in Region I and the coexistence line. This suggests that DWT gives a 



correct description of the relaxational behaviour at late times throughout these regimes. 

In contrast, the eigenvalue of the transition matrix determined from DWT does not 
coincide with our results in Regions II-IV [[J This means that while the shock profile 
considered in [45] remains the exact stationary state along the curve abq = 1, the slowest 
relaxational mode is no longer given by the particular implementation of DWT proposed 
in [26,27]. An obvious question is whether it is possible to reproduce our findings by a 
suitably modified DWT. 

To this end it is useful to consider our results for the eigenvalue of the lowest excited 
state as a function of b for fixed a and q. A particular example (a = 1.5 and q = 0.6) is 
shown in Fig. |4j 



Si 

0.002 h- 



0.5 \_--' 1.0 



Figure 4. The lowest excitation for q = 0.6 and a = 1.5. The functions defined 



in (3.2), (3.91 (dotted) and fcM (dashed) are both displayed for < b < 2. The 



excitation gap is a combination of (3.2) and (3.6 1, and drawn as the bold curve (red 
online). 



When b > a we are in the high density phase and the system is gapped. The 



corresponding gap is finite and given by (3.2). Decreasing b we are approaching the 



§ It was shown in Ref. [45] that exact multi domain wall solutions of the master equation exist more 
generally along the curves ab = q~ n for n — 1, 2, 3, . . ., but their precise properties have been analyzed 
only for n = 1 [46] . 

|| For the case of totally asymmetric diffusion this was already observed numerically in Ref. [23] and 
analytically in Ref. [31]. 
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coexistence line b = a, where the gap vanishes and the relaxation is purely diffusive. 
Decreasing b further drives us into the low density phase and the gap is again finite. We 
expect the gap to grow as b decreases. However, at b c = a' 1 ^ 3 (~ 0.87 in the example 
shown in Fig. |4|), the slope of (3.2) as a function of b vanishes and for b < b c (3.2) 



increases with decreasing b. For values of b smaller than the crossover point 5 C , the gap 



is no longer described by (3.2) but by (3.6), and remains constant 



It is now straightforward to reproduce these results within the framework of an 



effective DWT. In Region I the DWT prediction [45] for the gap coincides with (3.2) 

(3.15) 



where D± are defined above. However, as we cross over into Region II the gap predicted 
by this DWT no longer agrees with the exact result. We therefore modify the DWT 
as follows. We postulate that in Region II the density p + ceases to depend on b and 
remains fixed at p^ s = b c /l + b c . Retaining the expressions (3.14) for the hopping rates 



of the domain wall one finds that the gap is then given by Si(p , p^ s ). The value of p^ s 
is determined from the requirement that 

d£ 1 (p-,p+) 



dp^ 



= 0. (3.16) 

By construction this modified DWT reproduces the exact result for the relaxation rate. 

The modification of the DWT becomes more involved if a > q~ 3 / 2 . In this case 
there are two crossovers as is shown in Fig. [5] for the particular example q = 0.8 and 
a = 3. The first crossover separating Regions I and III takes place at b cl = q 2 a 1.92 




Figure 5. Eigenvalue of the lowest excitation for q = 0.8 and a — 3. The functions 



defined in ^2$, <Q9j (dotted) and ( |3.12[ ) (dashed) are displayed for < b < 4. The 
excitation gap is a combination of all three, and drawn as the bold curve (red online) . 



for the example shown in Fig. |5J). For values of b < b Ct \ the gap is given by (3.9). In this 
case the 0{L~ 2 ) correction to S± vanishes. A second crossover, now between Regions 
III and IV, occurs at b r ? = q 1 ^ 2 



d£ 1 (p-,p+) 



0.89 in the example shown in Fig. |5J) where 
= 0. (3.17) 



=PcB 



For values of b < b C 2 the relaxation rate is given by (3.9) and no longer depends on b 
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It is clearly possible to reproduce the exact relaxation rates by adjusting the 
densities p in the DWT accordingly. Unlike above, in this case we do not have 
a convincing heuristic argument for the first crossover. The modifications of DWT 
described here are completely ad hoc, all we can say is that at the boundaries between 
the various regions, levels cross and the precise nature of the relaxational dynamics 
changes. It would be very interesting to investigate whether in Regions II- IV the relevant 
excited states are still domain walls. 



4. Analysis of the Bethe ansatz equations 



In the following we derive the results summarised in the previous section. To this end 



we analyse (|2.1|) and (|2.3|) in the limit or large lattice lengths L. It is convenient to 

qzf\ 

z + a 1 + c/V 



introduce functions 

9b(z) 



In 



In 



*(1 . 
z(l — q 2 z 2 



(4-1) 



1 



In 



1 + qaz 1 + qcz 



+ In 



z + b l + d/z 



1 + qbz 1 + qdz / 

The central object of our analysis is the "counting function" [39-41], 

1 1 

iY L {z) = g{z) + -g b (z) + - ^ K(z h z), 



where K(w, z) is given by 

K(w, z) — — In 



w — qz 1 — q 2 zw 
1 — qw/z 1 — wz 



(4.2) 



(4.3) 



(4.4) 



Using the counting function, the Bethe ansatz equations (2.3) can be cast in logarithmic 
form as 

2tt 

— Ij , j = l,...,L-l. 



y L {z 3 ) 



(4.5) 



Here Ij are integer numbers. Each set of integers {Ij\ j — 1, . . . , L— 1} in (4.5) specifies 



a particular (excited) eigenstate of the transition matrix. Based on numerical solutions 



of (4.5) using standard root finding techniques, we assume that the first excited state 



always corresponds to the same set of integers 

Ij = -L/2+j for j = !,...,£-!. 



(4.6) 



The corresponding roots lie on a simple curve in the complex plane, which approaches 
a closed contour as L —>■ oo. The latter fact is more easily appreciated by considering 
the locus of reciprocal roots zj 1 rather than the locus of roots Zj. In Fig. [6] we present 

2, c = K" = -0.01, d = Kg S = -0.023, q = 0.1 and 



results for a 
L = 150. 



a,7 



5, b 



13,5 
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Figure 6. Root distribution and reciprocal root distributions for a — k+ 7 = 5, 
b = k+ =2,c= K~ = -0.01, d = Kg s = -0.023, q = 0.1 and L = 150 



In order to compute the exact large L asymptotics of the spectral gap, we derive 
an integro-differential equation for the counting function Yl(z) in the limit L — > oo. As 
a simple consequence of the residue theorem we can write 

T E fa) = <f lL f( z ) Y L{*) cot (\lY l {z)\ , (4.7) 

where C = C\ + C2 is a contour enclosing all the roots Zj, C\ being the "interior" and 
C*2 the "exterior" part, see Fig. [7| The contours C\ and C2 intersect in appropriately 
chosen points ^ and £*. It is convenient to fix the end points ^ and ^* by the requirement 

liCO = + 1l(0=vt-^. (4.8) 

Using ( 4.8[ ) in (4.3) we obtain a nonlinear integro-differential equation for the counting 
function Yl(z). Our goal is to solve this equation for large lattice lengths L through 
an expansion in inverse powers of L. In order to do so we first rewrite eq:sum2int by 
separating the contributions coming from C\ and C<i- When doing this it is useful to 
note that on the contour of integration we have by definition of the counting function 
that \mY L {z) = 0. As a result the imaginary part of Y L {z) is positive on C\ and negative 
on Ci- Using the fact that integration from £* to £ over the contour formed by the roots 
is equal to half that over C2 — C\ we find, 



iY L (z)=g(z) + ±-g h (z) + ^ / Kltr.:)Y^tr)dtr 



2tt J Ci 1 - e" lLy i(™) 2vr J C2 e lLY ^ w ) 

where we have chosen the branch cut of K (w, z) to lie along the negative real axis. 

Our strategy is to solve the integro-differential equation (4.9) by iteration. Once we 
have constructed the desired solution Y^{z) we determine the corresponding eigenvalue 
of the transition matrix from equation (2.1) by turning the sum over roots into an 
integral using (4.7) 

£= -£0-7^ / e(z)Y' T .(z)dz 

Z7T 

_L_t s(z)YHz) Lf 4z)YL(z) 
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Figure 7. Sketch of the contour of integration C in (4.7 1. The open dots correspond 



to the roots Zj and £ is chosen close to zl-i and avoiding poles of cot(L>L(z)/2). 



Here the constant Sq was denned previously in (2.8) and the "bare energy" s(z) is 

(q-l) 2 z .. , ( 1 1 



eiz 



[l-z)(qz-l) 



q) 



z — 1 qz — 1 



(4.11) 



5. Low and High Density Phases 

The low density phase for q < 1 is characterized by a > 1 and b < 1, while the high 
density phase corresponds to a < 1 and 6 > 1. In these phases we find that the locations 
of the end points £ and £* are such that a straightforward expansion of the correction 



term in (4.9) in inverse powers of L is possible (see e.g. [42,43] and Appendix B). The 
result is 



lYriz] 



9W 



1 If* 

g h (z) + — I K{w,z)Y[{w)dw 



L 



2tt 



tt (K'{?,z) K'(£,z) 



+ 0{L~ 



(5.1) 



12L2 V Yl(?) y l(0 
where the derivatives of K are with respect to the first argument. We note that here we 
have implicitly assumed that Y' L (C) is nonzero and of order O(L ). The integral from 
to £ is along the contour formed by the roots, see Fig. [7| In order to utilize complex 
analysis techniques it is useful to extend the integration contour beyond the endpoints, 
so that it pinches the negative real axis at points zf = z c ± iO (z c G E). This leads to 
the following expression 



l Y L {z)=g{z) + jg h {z) + ±- 



K(w, z)Y' L {w)dw 



— i K(w,z)Y[(w)dw + — [ K{w,z)Y[{w)dw 
2tt J v 2tt J z + 



YL(0 



(5.2) 



The key to the solution of (5.2) is that all terms have simple expansions in inverse 
powers of L. In order to find the eigenvalue (4.10) up to order 0(L~ 2 ) we need to solve 
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gg| to order 0(L- 3 ) 
Yl(z) = 



Substituting the expansions 

oo oo 



(5.3) 



n=0 



n=l 



back into (5.2) yields a hierarchy of integro-differential equations for the functions y n (z) 

. + 

1 



Vn(z) = 9n(z) + 



2vri 



K(w, z)y' n (w) dw. 



(5.4) 



The integral is along the closed contour following the locus of the roots, see Fig. [7j The 
first few driving terms g n {z) are given by 



9o(z)= -ig(z), 

9i{z) = -igb(z) + Ki + \iK(z c , z), 

g 2 (z) = k 2 + X 2 K(z c , z) + fx 2 K'(z c , z), 

g 3 (z) = k 3 + \ 3 K(z c , z) + fx 3 K'(z c , z) + u 3 K"(z c , z). 



(5.5) 



The functions g and g^ are defined in (4.1) and (4.2) and 

1 — qz c z~ l 1 — 



K(z c , z) = - ln(-z c ) + In 



ZZ r 



q 2 zz c 



(5.6) 



1 — qzz~ l 1 

The terms involving the kernel and its derivatives arise from Taylor-expanding the 
integrands in the integrals from to z~ and from z^ to £. Concomitantly the coefficients 
K n , A n , fi n and v n are given in terms of S n , r) n defined by (5.3), and by derivatives of 
y n evaluated at z c . Explicit expressions are presented in Appendix C We show how 
to construct a general solution of the set of equations (5.4) for n < 3 under certain 
restrictions on the values of the parameter a, b, c and d in Appendix A Having this 
solution in hand, we may determine the coefficients K n , X n , \i n and v n as follows. 
Substituting the expansions (5.3) into the boundary condition (4.8), which fixes the 
endpoints £ and £*, we obtain a hierarchy of conditions for y n (z c ), e.g. 

y l (0 =vo{t) + \yi{0 + ^MO + -.. 



2/o (2 C ) + t iyi( z c) + y'o( z c)($i + m)] + ■■■ 



L 



71 



TV — 



(5.7) 



Solving this equation order by order, we find that in all cases considered in the present 
work 



A 3 

^3 



fJ-2 = 
Z c (^3 



X 2 = Ki 

Ai 



0. 



Z c (z c — 1) K 2 



7r 2 X 1 (l + X 2 1 )z 2 c (l 



(5. 



2i"" v '" ~' 6 (1 + ^c) 2 ' 

Here the parameter Ai depends on the values of a and b. As will be described in detail 
below, it determines various crossover regimes within the low and high density phase. 
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Having determined the counting function we may use equation (4.10) to evaluate 



the corresponding eigenvalue of the transition matrix. Evaluating the necessary integrals 
in the same way as for the counting function itself we obtain 

L 
2^ 



£ 



I E(z)Y[(z)dz-lJ2 e n L ' n 

J Zc n>0 



(5.9) 



where the integral is over the closed contour on which the roots lie, e is given in (4.11) 
and 



e = \\e(z c ), 

e x = \2e(z c ) + n 2 s'(z c ), 

e 2 = \3e(z c ) + Hze\z c ) + v 3 e"{z c ). 



(5.10) 



Substituting the expansion for Yl(z) in inverse powers of L into (5.9) we arrive at the 
following result for the eigenvalue of the transition matrix with the largest non-zero real 
part 

I - 1A1 



£i=-(l-q) 



+ 



1 



2Zr 



1 + a ' 1 + 6 

>2W2 



1 



El- 



+ i (i- g )iAi(i + A;v + 0{L _ 3] 



L 2 



6(z„ 1 — z. 



(5.11) 



Here, the sum over m is over all poles of g' h {z) other than 0, — c and —d, that lie inside 
the contour of integration. The constants c m are the corresponding residues. We note 
that the number and position of such poles depend on the values of the parameters a 
and b. The values of both Ai and z c in turn depend on these poles. In particular we 
find 



Ai = 2i + i J^. 



(5.12) 



The result (5.11 ) for the smallest relaxation rate is generically a constant of order O(L ) 



implying an exponentially fast relaxation to the stationary state at large times. We note 



that due to the symmetry of the root distribution corresponding to (4.6) under complex 



conjugation, £\ is in fact real, and hence there are no oscillations in the slowest relaxation 
mode. 



5.1. Region I: large values of a and b 

The first regime we consider is obtained loosely speaking by taking q to be small, a and 
b large and positive, c and d small and negative. More precisely we require 



(i) —a and —b lie outside the contour of integration in (5.2) 

(ii) —1/qa and —l/qb lie outside the contour, 

(iii) — c and —d lie inside the contour of integration. 
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We will assume that the last assumption is fulfilled, postponing a detailed discussion 
to Section [6] Condition (i) amounts to the inequalities —a < z c and —b < z c , which 
translate to 

b > &d = a~ 1/3 (a > 1), a > o Cj i = r 1/3 (& > 1). (5.13) 

Condition (ii) implies that —1/qa < z c and —1/qb < z c , resulting in 

b<b c , 2 = q 2 a (a > 1), a < a c>2 = q 2 b {b>l). (5.14) 

From the distribution of the reciprocal roots, see Fig. |6j we infer that for these values 
of the parameters, the roots lie in fact inside the unit circle. We therefore assume, 
and verify a posteriori, that z c ^ — 1 and that the points ±1 lie outside the contour of 
integration. Combining the above assumption we conclude that the driving term (5.5) 
1 can be represented in the form 

z + cf 



for (5.4) with n 

9i(z) 



l In z — l In 



z + c 



lln 



+ Ailn 



z — qz c 



(5.15) 



where gf(z) is analytic inside the contour of integration. Under the above assumptions 
we may now solve the system (5.4) of integro-differential equations and then verify a 
posteriori that all underlying assumptions in fact hold. Some details of this caculation 
are presented in Appendix A The result for the counting function is Yl(z) = yo(z) + 

iyi{z) + t?V2{z) 



yo{z) 



= — i In 

= — l In 

— l In 

— l In 
+ Ailn 



0(L" 4 ) where 



z 

Zc 



1 - z 



(5.16) 



z 1 - 


z 2 ' 

^c 


+ Ki — i In 


(ab) — Ai ln(— z c ) 


Z c 1 - 


Z 2 _ 


(-c/z; 


q)oc 


(-cz; q)oo(- 


~z/a; q) OQ (-qaz c ; q)^ 


{-c/z c \ 


?)oc 


(-cz c ;q)oo{ 


-z c /a; q)oo(-qaz;q) 00 


(-d/z; 


q)oo{-dz;q) 00 (- 


-z/b] q)oo(-qbz c ; q)oo 



-d/z c ; q)oo(-dz c ; q) 
{qzc/z; q) 00 (qzz c ] q) 



y 2 (z) = k 2 - A 2 ln(-2; c ) - 



(qz/z c ;q)oo(qz 2 ;q) 

1'2 



■z c /b; q) oc ,(-qbz; q) 

oo 

1 

+ Ailn 



(5.17) 



+ fJ>2 

+ A 2 ln 

, ^2 



Mz\q) + 2ij 1 (z- 1 \q- l )\ 



i\z\q 



(qz c /z;q) 00 (qzz c ;q) 



2 1 

CO 



(qz/z c ;q)oo(qz 2 ;q) 
1 1 



+ A 2 ln 



z — z„ 
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+ A 3 ln 
+ ' 



{qz/z c ;q) 00 {qzl;q) 

fJ>3 ^3 A ( 1 



)oo 

oo 



+ A 3 ln 



z - z. 



-i 



Zr Z r . 



C C / \ C C C 



»3 



(z - z- 1 ) (z c - z c i y 



Here (a; q)^ denotes the q-Pochhammer symbol 

oo 

(a;q)oo = Y[(l-aq k ), 

k=0 

and we have defined functions 

oo 



n=0 



z c — q z) k z*{l — q 



-,n+l\k ' 



(5.19) 



(5.20) 



(5.21) 



Imposing the boundary conditions ( 4.8[ ), (5.7) and using the expressions presented in 

(5.22) 



Appendix C for the various constants we obtain 

1 



Ai = 2i, 



ab 



and 



A3 = A*2 = A 2 = «i = 0, 

/ \ 2 2 ( ) 
^3 = ^3 = ^c(^c - 1) ft 2 = -ITT -7— r^-. 



(5.23) 



This is in agreement with our previous assertion (5.8). Given our result for the counting 



function we may then determine the corresponding eigenvalue of the transition matrix 

2z r 



from (5.9) 



£1 = - a - ?) ( -!- + 1 



+ 



1 + a 1 + b 1- z c 
1 (l-g)Tr 2 



(zr 1 - 



+ 0(L~ 



(5.24) 



which is the result given in (3.2) 



5.2. Region II: —b inside the contour 

We now consider the case where —b moves inside the contour, but where the pole at 
—1/qa remains outside, i.e. 



b < b r 



-1/3 



(o>l). 



(5.25) 



The case where a < a c is readily obtained from the results below by the interchange 
a «-> b. The main difference compared to Region I is that the driving term g\{z) acquires 
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an additional branch point inside the contour of integration. The 0(L _1 ) contribution 
yi(z) to the counting function must therefore be determined on the basis of a different 
analytic structure of the driving term gi(z): 



gi{z) — — \ \nz — lln 
+ Ailn 



z + c 



— lln 



z + d 



qz c 



-iln(z + 6) +gl{z) 



(5.26) 



where g\{z) is analytic inside the contour. The solution of the corresponding integro- 
differential equation proceeds along the same lines as before, resulting in the expression 



(A. 29) for yi(z). The solutions of the equations for y 2 (z) and 2/3(2) remain unchanged. 
Imposing the boundary conditions (4.8) imposes Ai = 3i and z c = — a _1//3 , resulting in 



the eigenvalue (3.6) 



5.3. Region IIP —1/qa inside the contour 

The next case we consider is when —1/qa lies inside and —b outside the integration 
contour, which occurs in the parameter regime q 1 ^ 2 < b < q 2 a. The driving term gi(z) 
of the integro-differential equation for y z (z) is expressed as 

' z + c\ , / z + cf 



l In 



+ Ailn 



lln 



z 

z - qz c 

z 



+ lln 



z + 



l/qa\ 



(5.27) 



where gf(z) is analytic inside the contour of integration. Proceeding as before we arrive 



at the result for yi(z) given in (A.34). The results for 1/2(2) and 2/3(2) are the same as 



before and are given in (5.18) and (5.19). Imposing the boundary conditions (4.8) fixes 
\\ = 1 and z c = —q/b, leading to the eigenvalue given in (3.9). 



5.4. Region IV: —1/qa and —b inside the contour 

The last case we consider is when both —1/qa and —b lie inside the contour of integration. 
This occurs when a > q~ 3 / 2 and b < q 1 ! 2 . We may express g\{z) in the form 

' z + cf 



1 In 



+ Ailn 



z + c 

z 
z 



lln 



iln(z + b) 



qz c 



+ lln 



z + 



l/qa\ 



(5.28) 



z J \ z J 

where gf(z) is again analytic inside the contour and then proceed as in the other cases. 



Solving the integro-differential equation for y%{z) results in (A.39). The results for y 2 (z 



and 2/3(2) are again given by (5.18) and (5.19) respectively. Imposing the boundary 



conditions (4.8) now gives Ai = 2i and z c = —q 1 / 2 and leads to the eigenvalue given in 
d3l2l. 
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6. Dependence on c and d 

In all calculations described above we have assumed that both — c and —d lie inside the 
contour of integration. This is the case if 

— c < z* and —d<z*, (6.1) 

where z* is the point where the contour crosses the positive real axis, i.e. the solution 
of the equation Yl(z*) = 0. In leading order this is determined by the solution of 

Z* (1 - Zc) 2 , v 

z c (1 - z*) 2 

This condition is easily solved for z* as a function of z c and using the explicit expressions 
for z c in Regions I-IV we obtain corresponding restrictions on the allowed values of c 
and d. For example, in Region I we find 



1 + iVab + ab - \ (1 + iVab + ab) 2 - Aab 

-d < z*{a, b) = . (6.3) 

2yab 



However, as we will now show, the results in the high and low density phases presented 



above have a somewhat larger realm of validity than suggested by (6.3). To that end 
let us consider the situation where —d is still inside the contour of integration, but — c 
is slightly larger than z*(a,b). Then the root distribution corresponding to the largest 
eigenvalue has the same set of integers as before 

Y L (z j ) = -n+ 2 ^-, (6.4) 

but the root distribution now has an isolated root ( lying outside the contour of 
integration. The position of the isolated root is 

C = -c + 0(e-" L ) i (6.5) 

where v > 0. In Fig. [8]the root distribution is depicted for a specific case with c = —0.2, 
and it can be seen that ( ~ — c. The value of z*(a, b) in this case is approximately 0.136. 



Figure 8. Root distributions for a = k+ 7 = 5, b = nt s 



2, c = K~ 



-0.2, 



d = k 



0,5 



-0.023, q = 0.1 and L = 50 



In order to turn the summation over this distribution of roots into an integral, the 
contour C 2 of Figure [7] has to be extended to include the isolated root (. This is achieved 
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by adding a contour that runs from z* + iO to £ + iO and back from £ — iO to z* — iO. 
Importantly the extra contour does not encircle any other poles or branch points of the 
counting function. Hence the analysis of this case is exactly the same as before, and the 
eigenvalue is again given by (5.24) 

The situation becomes more complicated when — c is increased further. Now several 
isolated roots may lie on the positive real axis outside the contour, close to points — q n c 
for nGff. Naively the same argument as for the case of a single isolated roots applies, 
but a detailed analysis of such cases is beyond the scope of this publication. 



7. Conclusions 



In this work we have analyzed the Bethe ansatz equations of the partially asymmetric 
exclusion process with open boundaries. We have focussed on the parameter regime 
corresponding to low and high density phases in the stationary state. We have 
determined the eigenvalue of the transition matrix with the largest non-zero real part, 
which characterizes the relaxation towards the stationary state at asymptotically late 
times. We found that both the low and high density phases are subdivided into several 
regimes, which are characterized by different relaxational behaviours. In the vicinity 
of the coexistence line which separates low and high density phases the relaxational 
behaviour can be understood in terms of diffusion of domain walls. In the other regimes 
such interpretations are still possible, but are not conclusively supported by the available 
results. A number of open questions remain. We have not studied the parameter regime 
corresponding to the maximal current phase at stationarity. Furthermore, our analysis 
has been restricted to the "forward bias regime", in which the injection/extraction 
rates at the boundaries are compatible with the bias present in the bulk. It would be 
interesting to extend our analysis to the maximum current phase as well as the "reverse 
bias regime". Another open issue is the precise physical nature of the relaxational 
mechanism sufficiently far away from the coexistence line. Is the relaxation still driven 
by diffusion of some kind of domain walls, or does another mechanism take over? This 
issue is particularly relevant in Region III, which occupies a large part of the phase 
diagram when q — > 1. For symmetric diffusion q — 1 the relaxation is known to 
be diffusive and the spectral gap scales as 0(L~ 2 ). Further interesting questions are 
whether the Bethe ansatz solution of the open XXZ chain can be used to calculate 
current fluctuations [47] and whether it is possible to determine correlation functions by 
means of Bethe ansatz [48]. Perhaps recent algebraic insights may offer new tools for 
further analysis. Such insights include the PASEP being an exceptional representation 
of the two-boundary Temperley-Lieb algebra [49] through the connection with the open 
spin chain [31], or the connection with tridiagonal algebras [50,51]. 
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Appendix A. Details on the calculation of the counting function 

In this Appendix we describe how to solve the integro-differential equation for the 
counting function. We start by deriving a simple identity that proves to be very useful 
in the subsequent analysis. Let C be the contour of integration from z~ to z+ defined 
by the locus of roots, c.f. Fig|7j Let D denote the interior of C and f'(z) be an analytic 
function in D. Then elementary considerations show that 



r \n(w - z) + f(w)) dw 

2tti J - \w - a J 



z — a 



In ( ) + f(z c ) - f(z) if a £ D, z G D, 



(A.1) 



ln(a - z c ) + f(z c ) - f(z) if a, z G D. 



Here, for every w, we have placed the branch cut of \n(w — z) along the line from z to 
z c so that the left hand side is a well defined contour integral. The identity eqdogint is 



useful for constructing solutions of the integro-differential equations (5.4). In particular, 
if y(z) is given by 

y(z) =y a (z) + J2 A n l Hz-z m ), (A.2) 

rn 

where dy a (z)/dz is analytic in D, then, for z close to the locus of the roots, 

^ j> K(w, z)y'(w)dw = y(qz) - y(z c ) + ^2 A m K(z m , z). (A.3) 



Here the kernel K(w,z) is given in (4.4) whose branch cut we take as above. The 
condition on z is that 1/z, l/q 2 z and z/q all lie outside the contour of integration [zjq 
lies outside the contour as we take z to be very close to the latter). 



Appendix A.l. Region I: large values of a and b 
We first consider the case 

(i) — c and —d are inside D, 

(ii) ±1, —a, —b, —1/qa and —1/qb are outside D. 

Loosely speaking this corresponds to small q, large positive a and b and small negative 
c and d. We further assume (and verify a posteriori) that D lies inside the unit circle, 
and hence that the point 1/qc and 1/qd also lie outside D. 
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We wish to solve the system of integro-differential equations (5.4): 

. + 

1 



Vn{Z) 



9n{Z) + 



2vri 



K {w,z)y' n {w) dw, 



(A.4) 



where the driving terms g n (z) are given in (5.5). We note that the equations are coupled 
as the driving terms of the equations for larger n depend on the solutions for smaller 



values of n. In the following we construct a solution of (A.4) and verify a posteriori that 
the above assumptions hold. 

Appendix A. 1.1. Equation for yo{z): 



The driving term of the leading equation is given by 

1 — qz 



9o(z) = -ig(z) = -lln 2 - 2iln 



1 - z 



-ilnz + g%{z). 



(A.5) 



Here g§{z) is analytic in D. We assume that yo(z) has the same analytic structure in 
D, i.e. 

y (z) = -ilnz + (A.6) 
where y^z) is analytic. Subsituting this ansatz into the integro-differential equation 



(A.4) for n = we obtain from eqdogintK, 

y%(z) = g%{z) - iK(0,z) + y {qz) - y {z c ). 



(A.7) 



Combining equations (A.6) and ( A.7[ ) we obtain a functional equation for y^(z) 

y*(z) - y a (qz) = -2iln (j^fj + i H~*o) ~ ( A - 8 ) 

The constant yo(z c ) is readily determined by setting z = 0. The resulting functional 
equation is then readily solved, giving 



y%(z) = 2iln(l-z) -lln 



:i - z c ? 



(A.9) 



The zeroeth order term in the expansion of the counting function is then found to be 



y (z) = -lln 



z ( 1 — z t 



z,\l-z 



(A.10) 



Appendix A. 1.2. Equation for yi(z) 



With the assumptions on q,a,b,c and d stated above the driving term (5.5) of the 



integro-differential equation (A.4) for n — 1 can be cast in the form 
gi (z) — — l In z — i In ( — 1 — l In 



+ Ax In 



z — qz c 



+ 9t(z) 



(A.11) 
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where gf{z) is analytic in D. The singularities of yi(z) inside the domain D can be 
inferred by trying to solve the integro-differential equation by iteration. It is then quickly 
seen that a branch point in yi(z) at z — —c produces a branch point at z = —qc, which 
in turn leads to a branch point at z = — q 2 c etc. This suggests the following ansatz for 

yi(z) = - i\n[(-c/z;q) 0O (-d/z;q) 0C \ + Ai \n(qz c /z; q)^ 

-ilnz + y^z), (A.12) 

where yf(z) is analytic inside D, and the q-Pochhammer symbol (a; g)oo was defined in 



fl5.2Qp . Substituting ( |A.12[ ) into <\AA} we find 

yl(z) = lln [(-c/qz] q) ao (-d/qz; q)^] - A x ln(z c /z; q)^ 
- lln [(1 + cz){\ + qcz){\ + dz){\ + qdz)] 
+ Ai In [(1 - qz c z)(l - q 2 z c z)] + g*{z) - y x {z c 
+ \\\i(-qz) +y l (qz). 

This leads to the following functional equation for yf(z) 

y{(z) = - lln [(1 + cz)(l + qcz){\ + dz)(l + qdz)] 
+ Ai In [(1 - qz c z)(l - q 2 z c z)] 
+ 9t(z) - yi(z c ) ±tt + yl(qz). 



(A.13) 



The value of yi(z c ) is easily determined by evaluating ( A. 14) at z = 

yi (z c ) = ±7T + ^(0). 
Using the fact that for two analytic functions u and v, the equation 

u(z) — u(qz) = v(z), 



(A.14) 
(A.15) 
(A.16) 



is solved by u(z) = ^2^ = Qv(q k z), it is now a straightforward matter to solve the 



functional equation (A.14), leading to the result given in equation (5.17). 



Appendix A. 1.3. Equation for yi (z) : 



The driving term (5.5) of the integro-differential equation (5.4) for n = 2 can be 



represented in the form 



, z — qz c \ q 
g 2 {z) = A 2 ln ( — ] -/i 2 - — +g%{z), 



(A.17) 

z J z — qz c 

where g^iz) is analytic in D. As in the case n = 1 the singularities of yi{z) inside D 
can be determined by attempting to solve the equation by iteration. This results in the 
ansatz 



q n+l 



y 2 (z) = A 2 \n{qz c / z; q)^ - — 



n=0 



z — q n+1 z c 



(A.18) 
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where y 2 (z) is analytic inside D. Substituting (A. 18) into (A. 4) we obtain a functional 
equation for y 2 (z) 



y$(z) - y 2 (qz) = k 2 - X 2 ki(-z c ) + A 2 In 

qz 



+ 



1 — qz c z 1 — z c z qz — z c 



[1 - z c z)(l - qz c z) 
1 - qz/ z c 



~ V2(Z C ). 



Evaluating (A. 19) at z = fixes the constant to be 
y 2 (z c ) = k 2 - A 2 ln(-z c ) - — . 



(A.19) 



(A.20) 



The functional equation (A.19) is then solved by elementary means, resulting in the 



expression for 2/2(2) given in equation (5.18). 



Appendix A. 1.4- Equation for y 3 (z): 



The driving term (5.5) of the integro-differential equation (5.4) for n = 3 can be 



represented in the form 

93(z) = A 3 In 



z — qz c 



/4r 



Q 



- "3- 



(A.21) 



z J z - qz c [z - qz c ) 

where g 3 (z) is analytic in D. Determining the singularities of y 3 {z) inside D by iterating 
the integro-differential equation now results in the ansatz 

q n+l 



y 3 (z) = A 3 \n(qz c / z\ q)^ -^3^ 



n=0 



z — q n+1 z c 



q 



2n+2 



n=0 



'z-q^z c y +yAZ) 



(A.22) 



where y 3 (z) is analytic inside D. Substituting (A.22) into (A. 4) for n = 3 gives a 
functional equation for y 3 (z) 

2/3O) - Vliqz) = k 3 - A 3 ln(-z c ) - y 3 (z c ) 

[1 - z c z)(l - qz c z) 



+ A 3 In 



1 - qz/z c 
qz z 



V3 



1 — qz c z 1 — z c z qz — z c 
2 



q 2 z 2 



z 



_(l-qz c z) 2 (1 - z c z) 2 (qz-z c ) 2 _ 



Evaluating (A. 23) at z = again fixed the constant 



y 3 (z c ) = k 3 - X 3 ln(-z c ) 



V3 + Yl 

Zc, 



(A.23) 



(A.24) 



Solving the functional equation (A.23) then results in the expression for y 3 (z) given in 



equation (5.19). 
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Appendix A. 2. Region II: —b inside the contour of integration 

The determination of y n (z) for n > 2 is exactly the same as in the previous section and 



in particular the expressions for y%{z) and y%{z) are unchanged and given by (5.18) and 



(5.19). 



Appendix A. 2.1. Equation for t/i(z) 



With —b lying inside D the driving term g\(z) may be expressed as 



— 2i In z — i In I — 1 — i In 



+ Ailn 



z — qz c 



l In 



z + b 



(A.25) 



z J \ z 

where gf u (z) is analytic in D. Proceeding as before, we arrive at the following ansatz 
for y x (z) 

2/1 00 = - i\n[(-c/z;q) 00 ( y -d/z;q) 00 ] + Ai ln(g2; c /z; g)^ 

- iln(-6/z; g)oo - 2x1ns + (A.26) 



where y\(z) is analytic in D. Substituting (A.26) into (A.4) we obtain the functional 
equation 



y*(z) = - iln [(1 + cz)(l + qcz){\ + dz)(l + qdz)] 

+ Ai In [(1 - qz c z)(l - q 2 z c z)] - y\{z c ) ± 2n 
- i In [(1 + bz){l + qbz)} + g\Az) + y?(gz). 



Evaluating (A.27) at z = again fixes the constant 



J/l(Zc) = ±27T + <7? n (0), 



(A.27) 
(A.28) 



and solving (A.27) then results in the following expression for yi(z) 

^2" 



yi{z) 



— 2iln 

— l In 

— l In 
+ Ailn 



z 

Z r 



i In 



+ Ki — iln (a) — Ai ln(— z c 



l-z 2 _ 

(-c/z; q) OQ (-cz; q) OQ (-z/a; q) 00 (-qaz c ; q)^ 
(-c/z c ; q)oo(-cz c ; q)oo(-z c /a; q) 00 (-qaz; q) c 
(-d/z; q) 00 (-dz; q) OQ (-b/z; q) 00 {-bz\ q)^ 
{-d/z c ; q)oa(-dz c ] q)ao(-b/z c ; q)oo{-bz c ] q) 



(qz c /z; q)oo(qzz c ; q) 



)oo . 



+ Ailn 



z — z r . 



Z r — Z, 



-1 



(A.29) 



. (qz/z c ;q) 00 (qz^;q)l 
Appendix A. 3. —1/qa inside the contour of integration 

The determination of y n {z) for n > 2 is exactly the same as before and in particular the 



expressions for 3/2(2) and y^(z) are unchanged and given by (5.18) and (5.19) 



Appendix A. 3.1. Equation for y\ (z) : 



Slowest relaxation mode of the PASEP with open boundaries 



26 



In this parameter regime gi(z) is expressed as 



1 In 



z + c 



z 



lln 



z + d 



z 



+ Alhl ,^ +lln (i±i^) +9 - IllW , 



(A.30) 



where glm{z) is analytic inside D. The appropriate ansatz for yx(z) takes the form 

Vi(z) = -i\n[(-c/z;q) 0O (-d/z;q) 0C \ + Ai \n(qz c /z; q)^ 

+ i\n(-l/qaz;q) 00 + yl(z), (A.31) 



where yf(z) is analytic inside D. Substituting (A.31) into (A. 4) we arrive at the 
functional equation 



y*(z) = - lln [(1 + cz)(l + qcz){\ + dz)(l + qdz)] 
+ Ai In [(1 - qz c z)(l - q 2 z c z)] - yi(z c ) 
+ i In [(1 + z/qa){\ + z/a)} + gl in (z) + y^qz). 



Evaluating (A. 32) at z = gives 

yi (O = <7i,ni (o), 

and we finally arrive at the following expression for yi(z) 



(A.32) 
(A.33) 



— i In 

— l In 

— l In 
+ A x ln 



+ Ki — lln ip/q) — Ai ln(— z c 



-c/z; q) 00 (-cz;q) 00 (-l/qaz c ;q) 00 (-z c /qa;q) c 
_{-c/z c ] q)oo(-cz c ; q) oc (-l/qaz; q) 00 (-z/qa; q) c 

(-d/z; q) 00 (-dz; q) 00 (-z/b; q)oo(-qbz c ; q)^ 
_(-d/z c ; q)oo{-dz c ; q)oo{-z c /b] q) oa (-qbz; q) 



{qz c /z; q) 00 (qzz c ; q) 



)oo 



+ Ailn 



z — z„ 



(A.34) 



(qz/z c ;q) 00 (qz^;q)l 
Appendix A. 4- —b and —1/qa inside the contour of integration 

The determination of y n {z) for n > 2 is exactly the same as before and in particular the 



expressions for 2/2(2) and yz(z) are unchanged and given by (5.18) and (5.19) 



Appendix A. 4-1- Equation for y%(z): 

Now the driving term is expressed as 
gi(z) = — lln z — iln' 



, , z + d\ . ( z + b 
1 In — 1 in 



+ Al ln( Z —^ ) +iln 



z J \ z 

z + l/qa\ 
z ) 



(A.35) 



Slowest relaxation mode of the PASEP with open boundaries 



27 



where g\ w (z) is analytic in D. Our ansatz for yi(z) then takes the form 

yi{z) = -ilnz-iki[(-c/z;q) 00 (-d/z;q) 00 ] + Ai \n(qz c /z; q)^ 

- 1 ln(-6/ z\ g)oo + 1 hz(-l/qaz; q)^ + yl(z), (A.36) 
where yf(z) is analytic. Substituting (A.36) into ( A.4[ ) we obtain the functional equation 
y{(z) = - lln [(1 + cz)(l + qcz){\ + dz)(l + qdz)} 
+ Ai In [(1 - qz c z)(l - q 2 z c z)] - y x (z c ) ± vr 
+ lln [(1 + z/qa)(l + z/a)] - lln [(1 + bz)(l + qbz)] 
+ 9iMz) + yi(qz). 
Fixing the constant as before gives 



Vl {z c ) = ±7r + g*(0), 



resulting in 



— l In 

— l In 

— l In 



z 

Z r . 





ri-*c 2 i 


— lln 






1 - z 2 _ 



+ K\ — lln (g) — Ai ln(— z c 



+ Ailn 



(-c/z; q) OQ (-cz; q) ^(-1 / 'qaz c ; q) 00 (-z c /qa; q) c 
_{-c/z c \ q) OQ {-cz c ; q) 00 (-l/qaz; q) co (-z/qa; q) c 

(-d/z; q) 00 (-dz; q) OQ (-b/z; q) 00 {-bz; q)^ 
_(-d/z c ; q)oo(-dz c ; q)oo{-b/z c ; q)oo(-bz c ; q) c 
(qz c /z;q) OQ (qzz c ;q) 
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Ailn 



(qz/z c ;q) 00 (qz^;q) 
Appendix B. Analysis of the Abel-Plana Formula 



(A.37) 
(A.38) 



(A.39) 



In this appendix we sketch how to extract the finite-size correction terms from the 
integral expression 

>V z u; 



/ N 1 / \ 1 

9{z) + i9h(z) + — 



+ 



1 



K(w,z)Yl(w) 1_ 
2 J r l-e-^H 2vr 



K(w, z)Y[(w) dw 

K{w,z)Yl{w) 



C 2 



,iLY L (w) 



dw. 



(B.l) 



The main contributions to the correction terms in (B.l ) comes from the vicinities of the 
endpoints £, £*. Along the contour C\ the imaginary part of the counting function is 
positive, lm(Y L (w)) > 0, whereas along the contour C 2 it is negative, lm(Y L (w)) < 0. 
As a result the integrands decay exponentially with respect to the distance from the 
endpoints. In the vicinity of £, we therefore expand 

Y L (w) = Y L (0 + YL(£)(w - + . . . . (B.2) 

Assuming that Y[(£) is 0(1) (an assumption that will be checked self-consistently), we 
find that the leading contribution for large L is given by 



1 

2^ 



Ci 



K{w,z)Yj{w) 

1 _ e -iLY L (w) 



dw 



mo 

2tt 



K(w, z) 



I + Q -iLY[ (0(to-0 



dw 



(B.3) 
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Here we have used the boundary conditions (4.8). Carrying out the analogous analysis 
for the integral along the contour C2, we arrive at the following expression for the leading 

K(w,z)Yl(w) 



contribution of the last two terms in (|B.l|) 
1 

i 



A 



K(w,z)Y[(w) 



Ci 



2nL 



e -iLY L (w) 
1 



dw + 



o 



1 + e a 



KU + 



1 
1 x 



C 2 



dw 



1 x 



dx 



(B.4) 



If the endpoints £, £* are such that we can Taylor-expand the kernels appearing in (B.4), 
we can simplify the expression further with the result 



7T K'foz 



X 







1 + e x 

- (e- 



dx - (£ - r 



r: 



12L2 Y[(0 ^ * r (R5) 
This is the leading Euler-MacLaurin correction term that occurs in the low and high 
density phases, see (5.2). The key in the above derivation was the ability to expand 

ix \ ix \ 

-LYl(0(z- l -0+ix 



KU 



LYm 



KU 



In 



LYl(0(z-i-0-ix_ (R6) 

in a power series in x. This is unproblematic as long as LY^(£)(z _1 — £) is large, which 
turns out to be the case in the low and high density phases as well as on the coexistence 
line. 



Appendix C. Expansion coefficients 



In this appendix we list the coefficients arising in the expansion (5.4), (5.5) of the integral 
equation for the counting function Yl(z). In the following list we abbreviate y' n {z c ) by 
y' n . We note that by definition 5 n and rj n are real quantities. 



Ai 
A 2 

«3 



71 



7T 



y'^h - y'A - -y'o(Sl - 



vl) 



71 



- y' $3 - y'i$2 - ytfi - 2/0(^2 - mm) - 2 y "^ ~ vl) 
-\fM%-H), 



(C.l) 
(C.2) 

(C.3) 

(C.4) 

(C.5) 

(C.6) 
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A 3 = - (» + y[v2 + y' 2 vi + yo(<W + 5 2 rii) + <Wy" 




(C.7) 



A*3 = - 2/0(^1^2 + 8 2 r]i) + yi5i?7i + y'M^ - 4) + 

7T O 





(C.8) 



(C.9) 
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